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ABSTRACT 


We present an algorithm to approximately solve certain stochastic nonlinear 
programs through sample-average approximations. The sample sizes in these 
approximations are selected by approximately solving optimal-control problems defined 
on a discrete-time dynamic system. The optimal-control problem seeks to minimize the 
computational effort required to reach a near-optimal objective value of the stochastic 
nonlinear program. Unknown control-problem parameters such as rate of convergence, 
computational effort per solver iteration, and optimal value of the program are estimated 
within a receding horizon framework as the algorithm progresses. The algorithm is 
illustrated with single-commodity and multi-commodity network flow models. Measured 
against the best alternative heuristic policy we consider for selecting sample sizes, the 
algorithm finds a near-optimal objective value on average up to 17% faster. Further, the 
optimal-control problem also leads to a 40% reduction in standard deviation of 
computing times over a set of independent runs of the algorithm on identical problem 
instances. When we modify the algorithm by selecting a policy heuristically in the first 
stage (only), we improve computing time, on average, by nearly 47% against the best 
heuristic policy considered, while reducing the standard deviation across the independent 
runs by 55%. 
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EXECUTIVE SUMMARY 


Optimization of stochastic programs is challenging in part because there is no 
elosed-form solution, and beeause solution algorithms tend to be eomputationally 
intensive. The diffieulty arises beeause the objeetive funetion and/or eonstraint 
funetions are given by expeetations that eannot be evaluated exaetly. Henee, 
approximations are usually employed in optimization algorithms to estimate sueh 
funetions. One elass of sueh algorithms uses sample-average approximations within a 
nonlinear approximating problem (AP). However, there is no straightforward means to 
determine a sample size for use in these algorithms. Most algorithms resort to using 
heuristic policies for determining an appropriate sample size. In addition to sample size, 
beeause the AP is nonlinear, a suitable number of iterations of a nonlinear programming 
solver must also be seleeted for solving it. Most often, the number of solver iterations is 
seleeted prior to ealeulations beginning. It is often diffieult in praetiee to seleet sample 
sizes and number of solver iterations that balanees aeeuraey and eomputational effort. 

We develop a diserete-time dynamie system, the optimal eontrol of whieh 
determines a poliey for sample size and a number of solver iterations for use in the AP. 
The optimal-eontrol problem seeks to minimize the eomputational effort required to 
reach a near-optimal objeetive value of the stoehastic nonlinear program. The optimal- 
eontrol problem is approximately solved within a reeeding horizon framework, allowing 
repeated estimation of unknown parameters. 

Empirieal studies are performed on nonlinear single and multiple-eommodity 
network-flow problems on a grid constructed of 50 nodes with 134 arcs. The optimal- 
control problem selects sample sizes and solver iterations that lead to near-optimal 
objective values in less time than alternative heuristie polieies in all instanees tested. 
Measured against the best alternative poliey we eonsider for seleeting sample sizes, the 
algorithm finds a near-optimal objeetive value on average up to 17% faster. Further, the 
optimal-eontrol problem also leads to a 40% reduetion in standard deviation of 
eomputing times over a set of independent runs of the algorithm on identieal problem 



instances. The unknown parameters in the optimal-eontrol problem may be poorly 
estimated prior to the first stage of the algorithm, which may result in a poor poliey for 
the first stage. When we modify the algorithm by seleeting a poliey heuristically in the 
first stage (only), we improve computing time, on average, by nearly 47% against the 
best heuristie poliey eonsidered, while redueing the standard deviation aeross the 
independent runs by more than one-half. 
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I. 


INTRODUCTION 


A. BACKGROUND 

Optimization of stochastic programs has become a focus of much research over 
the past deeade. These problems are of partieular interest in part because they have no 
closed-form solutions and solution algorithms tend to be computationally intensive. The 
diffieulty arises because the objective function and/or constraint functions are given by 
expeetations that cannot be evaluated exactly. Hence, approximations are usually 
employed to estimate such functions. Approximations may provide precise, high-fidelity 
estimates with high computational cost, or may provide low fidelity with less 
computational cost. This research focuses on finding a balance between these two 
important faetors. 

Numerous design and planning applications require the optimization of 
stoehastie-programming problems. In aerodynamies, various problems arise in the 
optimization of a three-dimensional wing design (Alexandrov et ah, 2001). In one civil- 
engineering discipline, problems arise from structural optimization of bridges or support 
structures subjeet to failure probability constraints (Polak and Royset, 2007; Royset and 
Polak, 2007). Such problems may seek to optimize the eross-sectional dimensions of a 
support eolumn consisting of a material of particular yield strength subjected to various 
bending moments. Structural loading of the support eolumn weighs heavily on design 
considerations due to inherent failure probabilities of the materials. Many additional 
examples of stochastic-programming problems can be found in the areas of logisties and 
supply-ehain planning. Numerous papers diseuss solution methodologies for optimal 
design of supply-chain management including, but not limited to productions lines, 
consumer demand, availability of raw materials, warehousing and transportation of 
goods. In Santoso et al. (2003), a proposed stochastic-programming model and solution 
algorithm are used to eompute high-quality solutions to large-seale stoehastic supply 
chain design problems. Poojari et al. (2006) presents a method to solve discrete resource 
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allocation problems in the presence of future uncertainties for supply-chain planning and 
Sox and Muckstadt (1997) describe a finite-horizon stoehastic optimization model for a 
stochastic lot-scheduling problem. 

One approach to approximately solve optimization problems defined in terms of 
expectations is to use Sample Average Approximations (SAA), e.g., (Ruszezynski and 
Shapiro, 2007). That is, one or more approximating problems (APs) are solved, 
problems that replace eaeh expectation with a standard sample average approximation. 
The AP may be viewed as a deterministie mathematieal program, a nonlinear 
mathematical program in our case. A very large sample size in the sample-average 
approximation will provide a preeise estimate of the expected-value function, but the 
computational effort required to solve the AP for this sample size may be far too high. 
As an alternative, this researeh examines the eomputational effort associated with solving 
a sequence of APs where an efficient sample size is ehosen at each stage of the sequence. 
Coarse approximations are made early and, as the ealeulations evolve, adaptive 
adjustments to the sample size are made, inereasing the precision of the results. This 
adaptive control strategy is intuitively appealing as gains towards optimality eome 
initially at low computational cost through coarse approximations, and fidelity is 
inereased as larger samples are used near the completion of the algorithm. 

In addition to finding an efficient sample size for use with ealeulations, we are 
also concerned with seleeting an efficient number of iterations for the chosen nonlinear 
programming (NLP) solver used in solving various instance of the AP. The number of 
solver iterations has a profound impact on computational cost as well. Our approach will 
determine an effieient number of iterations for the NLP solver to carry out along with an 
efficient sample sizes for defining the APs as the overall algorithm progresses. 

B. LITERATURE REVIEW 

Stoehastic programming solution methodologies ineorporate both mathematieal- 
programming techniques and statistieally motivated approximations. These 
approximations are generally internal or external in nature. The distinction is found in 
the management of sample selection. With external methodologies (Royset and Polak, 
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2007), samples are taken before the solution algorithm begins (or they eould be taken 
before the algorithm begins), and no additional sampling is performed during the 
optimization proeess. In internal methods (Higle and Zhao, 2004), samples are intrinsie 
to the iterative solution proeess, performed whenever the algorithm requires the 
estimation of expeetations. An alternative to an external sampling approaeh with a fixed 
sample size is to vary the sample size during the ealeulations using elosed-loop or open- 
loop teehniques (Polak and Royset, 2007). With a elosed-loop teehnique, the sample size 
is adapted during the ealeulations. For example, if the objeetive value in a given iteration 
falls below a moving floor, the sample size may be inereased. On the other hand, in an 
open-loop teehnique, sample sizes are preassigned (Polak and Royset, 2007). In most 
oases, the sample sizes inorease as iterations progress towards an optimal solution. 

He and Polak (1990) desoribe a method for handling progressively finer stages of 
disoretization for semi-infinite optimization problems, whioh uses a preoision-adjustment 
strategy relevant to the present work. They formulate an auxiliary optimization problem 
that determines the number of solver iterations and preoision level at different stages of 
the ealeulations so that overall eomputing time is minimized approximately. The present 
work is motivated by this study. 

C. RESEARCH GOAL 

The goal of this researeh is to find an effieient way of seleeting sample size and 
number of solver iterations for use in solving a nonlinear stoehastie program through the 
solution of a sequenee of APs that use sample average approximations. We formulate a 
diserete-time optimal-eontrol problem that we solve approximately to obtain a preoision- 
adjustment polioy for determining the sample size and number of solver iterations. This 
polioy makes ooarse approximations in the early stages of the problem and, as 
progression to the optimal objeetive value is aohieved, the sample size inoreases to ensure 
oonvergenoe to a looally optimal solution. The optimal-eontrol problem is formulated 
with the objeetive to minimize the amount of eomputational work neeessary to reaeh a 
looally optimal solution to the AP. 
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D. STRUCTURE OF THESIS AND CHAPTER OUTLINE 


This thesis is organized into five chapters including the Introduction. Chapter II 
provides a discussion on the formulation of a nonlinear stochastic problem, an outline of 
the algorithm we intend to use to solve stochastic nonlinear network flow problems and 
introduces our precision-adjustment problem. Chapter III introduces the methodology we 
propose to use for selecting an efficient sample size and number of solver iterations in the 
solution algorithm; stopping criterion are also discussed. Chapter IV formulates and 
solves two types of stochastic nonlinear network flow problems. It compares our 
algorithm’s efficiency with the efficiency of similar algorithms that use heuristic sample- 
size and solver-iteration policies. Chapter V summarizes the research and presents the 
main findings and insights. 
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II. PROBLEM DEFINITION AND FORMULATION 


This chapter formulates a stoehastic nonlinear program and presents a eoneeptual 
algorithm we use as our solution approaeh. Further, it deseribes the preeision-adjustment 
problem we faee within the eoneeptual algorithm. 

A. PROBLEM FORMULATION 

This seetion deseribes the formulation of the “expected-value problem” and an 
approximating problem whieh will be used as a surrogate to approximately solve the 
expeeted-value problem. 

1. Expected-value Problem 

Consider the optimization of, say, an engineering design or logisties problem 
defined by the random funetion F{x, co ), where x e A c j " is a veetor of eontinuous 
deeision variables and ni is a veetor of eontinuous random variables defined on the 
probability spaee (Q, F, P). This situation results in a diffieult stoehastie optimization 
problem of the form; 

EP : min {/(x) := E[F(x, co )]}, (2.1) 

where X is a eompaet subset of \ " representing feasible solutions of the problem, and E 
is the expeetation taken with respeet to the known probability distribution P of the 
random veetor co. We assume that, for every x e X , the expeeted-value funetion is well 
defined and smooth (eontinuously differentiable). Eurther, we denote the optimal value 
of (2.1) as f* with a set of optimal solutions denoted by X*. Eurther we indieate a set of 
£-optimal solutions by Xj, i.e., for any ^ > 0 

X;:={xEX|/(x)-r <^}. (2.2) 

Beeause the distribution of F{x,co) is unknown, we eannot eompute the 
expeetation in elosed form, so we approximate the expeeted-value funetion by the 
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sample-average estimator diseussed below. Uneertainty may also be introdueed in the 
eonstraint funetions defining X; however, our researeh foeuses only on uneertainty in the 
objeetive funetion. 

2. Approximating Problem 

To approximate /(x) = E[F(x, (u)] we use a sample average defined by: 

(2.3) 

where O) ,o}^,..., 0 /, is a sample of size A eonsisting of independent, identieally 
distributed (iid) random variables. Moreover, we define an approximating problem (AP): 

AP: min/^(x), (2.4) 

xeX 

with an optimal value denoted by /^. We refer to (2.1) and (2.4) as the “expeeted- 
value” and “approximating problems,” respeetively. 

The expeeted-value funetion, for instanee, may be represented by the funetion 
depieted by the solid line in Figure 1. This funetion, as stated above, eannot be eomputed 
in elosed form, so we approximate the funetion’s form by the sample average as shown in 
Figure 1. As discussed in Section 3 below, we expect the sample average to take the 
form of the expected-value function as the sample size approaches infinity. It is known 
that solving the AP with an appropriate sample size provides a reasonable approximation 
to the solution to FP (Ruszczynski and Shapiro, 2007). 



Figure 1. Expected Value function compared to Sample Average function 
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3. Properties of Sample Average 


Monte Carlo simulation can be used to generate an lid sample of (u's of size N to 
obtain F{x,co^),F{x,co^),...,F{x,co^) for use in (2.3). It is well known that /^(x) is an 
unbiased estimator of /(x) i.e., 


E[A(x)] = /(x), 

Var[/^(x)] = cr"(x)/A, 

where 

cr^ (x) = Var(F(x, co)) . 

Because the generated sample is iid and given rather weak general assumptions, it 
follows from the Law of Large Numbers (LLN) that /^(x) converges pointwise to /(x) 

with probability one, as A ^ oo (Ruszczynski and Shapiro, 2007). Therefore, we would 
expect the optimal value and optimal solution of the AP to converge to those of EP, as 
A ^ 00 . This is made precise in the following proposition. 

Proposition 1. (Prop. 4.2; Ruszczynski and Shapiro, 2007). If the pointwise 
LLN holds, i.e., /^(x) converges to /(x) uniformly on X, with probability one 

as N ^ CO , then f* converges to f* with probability one, as A ^ oo . 

Moreover, if E[F(x,(»)^]< oo , then it follows under weak assumption from the 
central limit theorem (CET) that Va(/^(x)-/(x)) ^ T^, where ^ denotes 
convergence in distribution and is a zero-mean normal random variable with variance 
cr^(x) (Ruszczynski and Shapiro, 2007). Hence, for a large A, /^(x) is approximately 
normally distributed with mean /(x), and variance cr^(x)/ A . Eigure 2 illustrates this 
result. 

We know that 


min 

XE.X 


E[/„(x)]>E 


min/„(i-) 

XE.X 
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and, because /^(x) is an unbiased estimator of f{x), it follows that < /*. That 

is, the AP’s optimal objective value is a downward-biased estimate of the optimal 
objective value of the EP. 



B. ALGORITHM 

We consider a conceptual algorithm of the following form for approximately 
solving EP. 

Conceptual Algorithm (CA) 

Data; Optimality tolerance e > 0; initial point Xg e X . 

Step 1: Set stage counter k=\. 

Step 2; Determine N^. and . 

Step 3; Carry out n^. iterations of a solver applied to; 

min A (x), started with Xg using the found at Step 2. 

Step 4; If x^^ e X* , then Stop. Else set x^^' = x^^, replace khy k+\ and 
go to Step 2. 
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We refer to Steps 2-4 as a “stage.” CA uses a single set of sample realizations for 
eaeti stage, but all samples are independent between stages. For purposes of this 
researeh, we have ehosen the Projeeted Gradient Method (PGM) as our nonlinear 
programming solver (for example, see Polak 1997, p. 66). However, it is worth 
mentioning that any linearly eonvergent nonlinear programming solver may be used. Our 
goal is to determine and for eaeh stage that approximately minimizes the total 
computational effort required by CA to reach a near-optimal objective value to EP. We 
refer to the rule for selecting N^. and as a “policy.” 

C. PRECISION-ADJUSTMENT PROBLEM 

We define the task of selecting the appropriate and n^. for each stage of the 

CA as the Precision-Adjustment Problem (PAP). We formulate PAP as a particular 
discrete-time optimal-control dynamic program, which is discussed in the next chapter. 
PAP is solved approximately using dynamic programming to obtain a precision- 
adjustment policy for determining A^ and . This policy adaptively adjusts the sample 
size and number of solver iterations between stages. PAP is formulated to minimize the 
amount of computational work necessary to reach a near-optimal objective value of the 
stochastic nonlinear program. After carrying out solver iterations, we abandon further 

progress to the locally optimal solution for the current AP, and use this iterate as a warm 
start for the next stage of CA. 
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III. METHODOLOGY 


This chapter describes how we control the precision within the conceptual 
algorithm and how we estimate the parameters for the precision adjustment. Further, the 
chapter presents the criterion by which we stop the conceptual algorithm (in Step 4). 

A. PRECISION ADJUSTMENT CONTROL 

This section presents a discrete-time dynamic system, the optimal control of 
which will determine a sample size and a number of solver iterations n^. to be used 
for each stage of CA. 

I. Controlling Sample Size 

The dynamic system described in Royset (2009) lays the foundation for this 
research. Following his approach, we first identify the asymptotic distributions of the 
progress made by CA. We assume that the solver used in Step 3 of CA is uniformly 
linearly convergent (see Royset, 2009), i.e., there exists a ^e(0,l) such that 

/^(R^(x))-/; <6’(/^(x)-/;) for all xeX and Ae¥, where ¥= (1,2,3,...}and 
Pf^{x) is the iterate found after carrying out one iteration of the solver used in the CA 
starting from x with sample size N. We refer to 6 as the rate-of-convergence 
coefficient. We let (x) denote the iterate from the solver found at stage k, after n^. 
iterations with a sample of size N^. starting from x. It follows from the assumption of 
linear eonvergence and from the optimality of , that for any x e X , 

/; ^ u (Pxi (^)) ^ /;+(A (^) - /;) 

with probability one. Moreover, based on rather weak assumptions and theorems 
introduced in Royset (2009), if andn^. are large, approximately 

distributed as 
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^(/* +(/(<:; )-n,c7\x) IN,), o. i) 

where JV{fj.,<j) represents a normally distributed random variable with mean jj, and 
varianee cr. Figure 3 illustrates a possible approximate distribution of as stated 

from this eonelusion. 



Further, Royset (2009) shows that we ean heuristieally approximate the 
distribution with truneation at /* to aeeount for the faet that 

/(x) > /* for all X e X . Therefore our approximate distribution for /(xf^ |)) 
eonditional on xf ' is: 

^k-\ 

mr+o"^ {f{xi 2 ^ )-r),a\x)/N„r), ( 3 . 2 ) 

where JV denotes a truneated normally distributed random variable based on 
X(/*+ ^"‘(/(xf^ J)-/*),cr^(x*)/X^), with a truneation threshold /*. Figure 4 

illustrates this distribution. As diseussed in Subseetion 2, we ean eontrol this distribution 
by ehanging N, and n,. 
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Figure 4. Approximate Truncated Normal Distribution for /) 

Since (3.2) approximately holds for any k, we can use that expression to estimate 
the quality of a solution generated by the CA after any number of stages given a selection 
of sample sizes and number of iterations. The situation is illustrated in Figure 5, where 
the /(x° ) and a selection of (Aj, n^) determines the probability distribution of 

Given an outcome of that distribution, a selection of (Aj, n^), gives the distribution of 

/(^i),etc. 

The values for /*, 6, and cr^(x*) in the distribution (3.2) are unknown, however. 
Since x*^ J is known at the beginning of the stage, we can estimate /(x^^ |), and we 
construct estimation schemes as shown in Subsection B below to 
estimate /*, 6, and cr^(x*). We now develop a dynamic program as discussed in Royset 
(2009). 
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/ {x\ ) is a truncated normal with distribution 

{rn,a\x)i N„n 


f (x^J is a truncated normal with distribution 

ir+e'’^inx\)-r),a\x)/N„r) 

f (x^ ) is a truncated normal with distribution 
_ {f+e''^{f{xi)-n,cT\x)i N„r) _ 

Figure 5. Truncated normal distributions of function values 

2. Dynamic Program 

Beginning at stage k, we define estimates of /*, 0, and cr^(x*) as 

r^, r*, fg, and r ^2 respectively and as in Royset (2009), use (3.2) as a basis for a model 
of l = k, k + l, k + 2,.... We let fi, l = k, k + l, k + 2,... denote our estimates of 

l = k,k + \,k + 2,.... Controls are defined as {N^,n,^), (A^^^i,n^^i), 
and the dynamic equation for the state f, is 

= J^(r*+r;'(/,-r*),r^2 !Nj,r ), I = k, k+ \, k + 2,... (3.3) 

with initial condition and where the equality indicates equality in distribution. 

We define c{N, n) to be the computational cost of carrying out n iterations of the 
solver applied to the AP with a sample of size N, where the terminal state is characterized 
by c(l,0) = 0. Terminal states discussed in Royset (2009) are defined by Aj, and are 
translated for our case as 

T;={^ei \^-r*<8}. 
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The set of feasible eontrols R{^) are defined as; If 
then = {(1,0)}; otherwise, i?(^) = ¥ x¥. We seek an admissible stationary 
poliey r: j ^¥ x¥ u{0} with T{f,)GR{f,) whieh minimizes the eomputational eost 
by evaluating a planning horizon of feasible eontrols , 

• The sample-size eontrol-problem stated in Royset (2009) whieh 
aeeomplishes this task is 




,rg,r^):= limsupE’ 


Sc(r(/;)) 


l=k 


subjeet to eonstraints (3.3). Here E is the expeetation with respeet to the disturbanees on 
the stages k, k + s due to the truneated normal distribution in (3.3) Finally, we 
define the surrogate sample-size control problem by 


S - SSCP, {rf,r\rg,rj: inf J,^^ir^,r\rg,r^). 



Figure 6. Representative illustration of dynamie program 


(3.4) 


Now, assume we have estimated f(Xg) = /(xf |) at the beginning of stage k; see 
Figure 6. We wish to ehoose a pair {N^,n^) that will be eomputationally effieient so that 
at the end of the k'* stage shown in Figure 5 we have progressed toward /* at a 

“reasonable” eomputational eost. As shown by the dotted line in Figure 5, we define a 
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tolerance s to determine a range of near-optimal objective values that are acceptable for 
selection as /*. Further, we assume that the computational cost of the stage is 
defined as c{Ni^,n^) = Nj^n^.. Then, from (3.2), we see that selections of and have 
varying effects. A small increases the variability in the truncated distribution 
whereas a large compresses the distribution, reducing the variability. The number of 
solver iterations also impacts the cost by affecting the mean of the truncated distribution 
in (3.2). As increases, the mean of the distribution moves closer to /*, but also 
increases the computational cost. The optimal policy of S-SSCP^(r^,r*,rg,r^) balances 
the computational cost of selecting large A^. and and the likelihood that the CA 

reaches a near-optimal objective value in the current stage. Hence, we expect that policy 
to be reasonably efficient even though it is based on several approximations. 

To solve S-SSCP^.(r^,r*,rg,r^) approximately, we discretize the state space and 

the truncated normal distribution and then apply backward recursion to the resulting 
dynamic program. We refer to the policy computed in this manner as Model-Predictive 
Control (MPC). We adopt the discretization technique, as in Royset (2009), on a 
planning horizon of 10 future stages. A general depiction of the discretization scheme 
applied to Figure 6 is shown in Figures 7 and 8. An example of a relatively low cost 
selection of A^ and n^. is shown in Figure 7. Here, a relatively small Nk with a choice of 

small Hk provides little gains toward f*. Although progress is made toward f *, 
additional stages of CA are necessary arrive at a near-optimal value. 
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Figure 8 depicts a situation where a large Nk and a large Uk provide a near-optimal 
objective value close to /* but at a high computational price. For this situation, future 
stages of CA may not be necessary in order to achieve a near-optimal objective value. 

B. PARAMETER ESTIMATION 

We rely on estimates of the parameters 0, and in 

S-SSCP,(r, ,r*,rg,r^) as we describe next. 
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1. Estimating Variance 


Upon completion of iterations of the solver with a sample of size N^. in stage 
k, the set of iterates and function values are known. We use this 

information to estimate parameters r^, r\ rg, and • First, we discuss our estimate for 
r ^2 • We let denote our estimate of cr^(x*) for stage k + \ and set it equal to the 
sample varianee at the last iterate, i.e.. 


r 

(7 


2 


= := 


N,- 


i 7=1 




2. Estimation of Rate-of-Convergence Coefficient and Optimal Objective 
Value 


Next, we determine , the estimate of the rate-of-convergenee eoefficient 0 of 

the nonlinear programming solver used in Step 3 of CA. We adopt the method in He and 
Polak (1990) to estimate 0, but modify it slightly, adding an exponential smoothing step 
to avoid large changes in the estimate. As the estimate of 0 is updated after each stage, 
we let 0/^ be the estimate of 0 available at the beginning of stage k and set in 

S-SSCP,(r,,U,r„0. 

As in Royset (2009), the (biased) estimate of /* at the beginning of the k‘^ stage 
is computed by a weighted average of estimates of / = 1,2,...,A:-1, and is denoted by 
//. Unlike the approach used by Royset, f* is computed with a fixed, conservative 

estimate of the rate of convergence denoted by 0 . The meaning of “conservative” is 
discussed in more detail in the explanation of Step 6 of Subroutine PE below. 

Subroutine PE(k. 0, f *, (xf )};*„) 

Parameters; Exponential smoothing parameter ^g( 0,1), conservative estimate 
0 G (0,1) of rate of convergence 0 and tolerance > 0 . 
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2 . 


Input data: Previous stage index k; estimates and //; function values 
{/ 7 v,(^f)}"=o from stage k. 

Output: Returns estimates and/^*^j. 

Step 1: Set 9 = 0^. 

Step 2: Estimate the minimum of 

1 




i=0 




Step 3: Solve least-square problem: 

% 

(a,b) = arg min ^ (log(/^^ (xf ) - h) - z log a - bf 


i=0 


Step 4: If I (9 - a |< set ^ = a and go to Step 5. Else, set 6 = a and go to Step 
Step 5: Set 6^^^ =\i/6 + {\- y/)6i^ 

Step 6: Compute conservative estimate of the minimum value of 


m, := mm 






\-9’ 


Step 7: Estimate /*: 




N,^ 


i:y> 


—7- m, -\ -—- C 


Step 8: Return 6^^^ 2 LnA f*^ 


Erom our assumption of a uniformly linearly convergent solver, it follows that 
/; >(A^(x;)-^”‘-X(xf))/(l-^"‘-') for all i = 0, 1, 2, ..., n,-\. Step 2 averages 
these lower-bounding estimates with the current estimate of 6 and uses this as an 
estimate of . With this estimate, we compute an estimate of the rate-of-convergence 
coefficient in a least-squares sense. That is, we use log-linear regression to compute the 
rate-of-convergence coefficient to best fit the sequence We define a 

tolerance for computing the rate-of-convergence coefficient, denoted by Sg , which 
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determines when to terminate the regression technique. Next, we perform exponential 
smoothing so as not to let the estimate of 0 vary too much from one stage to another. 
The value in Step 6 of Subroutine PE is guaranteed to be a lower bound on only 

if 6* > 6*. Hence, we recommend 6 be set to a value close to 1. Since is simply the 
weighted average of m,, / = 1,2,...,A:, and E[/^]</* (see Chapter II, Section A, 
Subsection 3), we therefore find that is a lower bound on/*, on average, when 

e>e. 


C. STOPPING CRITERION 

We could implement a stopping criterion based upon a hypothesis test of Karush- 
Kuhn-Tucker (KKT) conditions as developed in Ruszczynski and Shapiro (2007). As 
stated in Ruszczynski and Shapiro (2007), suppose that the feasible set X is defined by 
equality and inequality constraints in the form 

A:={x€ i " :g.(x) = 0, i = \,...,q\ g.(x)<0, 1 = ^ + 1,...;?}, 
where the g,(x) are smooth deterministic functions. If only equality constraints are 
present and the gradient vectors Vg.(x), i = \,...,q are linearly independent, then the 

hypothesis test of KKT conditions can be based upon an asymptotically noncentral chi- 
square distribution. If the assumption of linearly independent gradient vectors cannot be 
met, a degenerate solution is presented due to redundancy in the constraint functions. In 
the case of both equality and inequality constraints, a similar result is available (see 
Ruszczynski and Shapiro, 2007), which also relies on the linear independence of 
gradients of active constraints. Since such linear-independence assumptions may often 
fail in practical application, we have chosen to adopt an approach using optimality gaps 
as defined by Mak et al. (1999), and similarly by Ruszczynski and Shapiro (2007). We 
base our stopping criterion upon this approach with a small modification. 
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With the iterate found at the completion of stage k, we again estimate 
as before, with /^. ) using a new independent sample of size N*. Here we elect to 

use a large sample size to obtain an accurate approximation of /(x^^) . Calculation effort 
is not substantially increased by this procedure as we are not performing an optimization. 
From the central limit theorem, a probabilistic upper bound on /* is approximately 
normally distributed with mean /(xf^) and variance )/ N* for large N". 


The modification of the method described in Mak et al. (1999) occurs in our 
construction of a lower bound as described in Subroutine PE. While Mak et al. (1999) 
use the average of the optimal values of a set of APs, we construct a lower bound on /* 
by averaging lower bounds on optimal values of the APs for each stage. Our method 
tends to be more conservative as it is based upon an assumption of a rate of convergence. 
From Royset (2009), we see that our lower bound is approximately normally distributed 

with mean /* and variance cr^(x*)/A,, for large Aj, Aj,..., and • 

Hence the inequality 


Prob(/(xf^)</*+f)> O 



(3.5) 


holds approximately. We therefore stop the calculations when the right-hand side in (3.5) 
exceeds a selected confidence level S , typically 0.95 or 0.99. 
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IV. COMPUTATIONAL STUDY 


This chapter presents a eomputational study of two network-flow problems to test 
how well Model-Predietive Control reduees eomputation times eompared to alternative 
polieies for seleeting sample sizes. We define two network-flow problems with random 
eongestion and present results of Model-Predietive Control as eompared to heuristie 
eontrol of sample sizes. 

A. SINGLE-COMMODITY NETWORK FLOW 


To develop a single-eommodity flow problem (SCF) for testing, we consider the 
generie eongestion model for single-eommodity flows as deseribed by Ahuja et ah, 
(1993, p. 651), but modify it to inelude random eongestion. Here, the generie model has 
a nonlinear objeetive funetion of the form 


Z 


X,J 


min 

x<^x —X 

(i,j)eA ij Xy 


(4.1) 


where M^j is the nominal eapacity of are (/, j) and Xy is the flow of a single eommodity 

on are For review of eommodity flow and eongestion modeling we refer the reader 
to pages 109-124 in Hearn et. al. (2001), and to Mareotte and Nguyen (1998) and 
Bergendorf et al. (1997). We first present an SCF problem and then advanee to a multi- 
eommodity flow problem. Even though the SCF problem is a speeial ease of the 
multiple-eommodity flow problem, we present SCF first, due to the relative ease of 
explaining this simpler problem. 


We eonsider a graph G = (N, A), where N and A are sets of nodes and ares, 
respeetively. The speeifie graph eonsidered in this study is shown in Figure 9. This grid 
network ean flow eommodity left-to-right, north-to-south and south-to-north, but not 
right-to-left. The test-problem grid has 50 nodes and 134 ares. The start node, denoted 
by 5, is the supply souree and the terminal node, t, is the demand sink. Individual eosts 
assoeiated with the flow aeross an are are assigned as random numbers from a normal 
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distribution with parameters that will be speeified, and we define a eongestion parameter 
on eaeh are by generating a log-normally distributed random variable, with parameters 
that will be speeified. 



Figure 9. Transportation Grid Network 


We formulate an SCF problem as follows: 

Indices 

i,j nodes, i,jeN 

(i,j) arcs (i,j)e A 

Data 

M.j eapaeity of are {i,j). 

C.j unit eost to ship eommodity on are {i,j). 

mean of log-normal random variable denoting eongestion on are {i,j). 
<jy varianee of log-normal random variable denoting eongestion on are 

D demand of eommodity. 
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Random Variables 


cl).j congestion parameter on arc (/, j); this is a log-normal random variable, 

with mean and variance afj. 

Decision Variables 

Xy amount of commodity shipped on arc 

Mathematical Formulation 


min 




y - 


s-t. i; X..- X X, 

f.(iJ)GA 


-D if z = 5 
< 0 if i^N\{s,t} 

D if i = t 


0<Xy<My, 


V 


hj 


(4.2) 


We note that the expeetation in the objeetive funetion ean be eomputed by 
evaluating 1^41 one-dimensional integrals, and thus a simpler method for solving this 
model is available. However, this model serves as a simple example to illustrate our 
solution approach, which applies to more general situations. 

We assign 500 units of supply at node 5, with a corresponding 500 units of 
demand at node t. Are eapaeities are chosen as =100 for all (z,y). Based on 

preliminary numerical experiments, we find that 6 = 0.993 is sufficient to obtain lower 
bounds on in Step 6 of Subroutine PE. 


B. MULTI-COMMODITY NETWORK FLOW 

For testing, we also eonsider a eongestion problem for a multi-commodity flow 
problem in a transportation network (MCE), and reuse the grid network described in 
Subsection A. The formulation is as follows: 
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Indices 

i,j nodes,/, 7 e 

O', 7 ) arcs (i,j)e A 

p commodity, 7 » e { 1 , 2 ,...,P} 

Data 

M.j capacity of arc 

unit cost to ship commodity p on arc 

p.. mean of log-normal random variable denoting congestion on arc 

<jy variance of log-normal random variable denoting congestion on arc 

D’’ demand of commodity p. 

Random Variables 

congestion parameter on arc (/, 7); this is a log-normal random variable, 
with mean and variance . 

Variables 

xf. amount of commodity p shipped on arc 

Mathematical Formulation 


min 


E 


I 


_ P _ 


P 



s.t. 


I I 

J:(J,i)eA i:(iJ)eA 


-D” if/ = 5 

< 0 if / e N\{s,t), V p 
D” ifz= t 


(4.3) 


V ij 
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0<x|, 


V i,j,p 


For this problem, arc capacity is increased to M.j =150 for all (i,j) to allow for 

increased flow from additional commodities. We consider two commodities with 
supplies at 5 equal to 500 and 300 and demands at t equal to 500 and 300, respectively. 
In MCF, preliminary experimentation shows that 6 = 0.997 tends to provide a valid 
lower bound of in Step 6 of Subroutine PE and we adopt that value for 6 . 

C. COMPUTATIONAL STUDY 

For this computational study, we apply the parameters described next to both SCF 
and MCF. We use the PGM nonlinear-programming algorithm with Armijo step size 
rule; for example, see, Polak (1997, p. 67) and Bertsekas (1999, p. 31). The quadratic 
direction-finding problem in the PGM is solved using LSSOL (Gill et ah, 1986) as 
implemented in TOMLAB 7.0 (Holmstrom, 2008). We use parameters 
a = 0.5 and (3 = 0.8 in the Armijo step-size rule and in Subroutine PE use an exponential 

smoothing parameter ^ = 1 / 3 and tolerance = 0.0001. 

Eor stopping criterion, we draw a new independent sample of size N* = 10000 to 
evaluate/^. (x^^) and use a stopping confidence level of (5 = 0.95. We use 

(Ay, cr,j) = (3,4) for all {i,j) as parameters for the log-normal distributed random variable 
a>^J representing congestion. Arc costs Cy and Cy are generated from a normal 
distribution of random numbers with mean 80 and standard deviation 20. Additionally, 
we set the relative optimality tolerance to 0.01 for use in calculations, i.e., s = 0.01 
on stage k. 

Eor comparison studies, we consider two versions of Model-Predictive Control; 
MPCl and MPC2. In MPCI, Model-Predictive Control is applied to all stages of the 
conceptual algorithm, including the first stage. We find empirically that MPC might 
have poor control in the initial stage when the parameters estimated in 
S-SSCPj.(A,r*,rg,r^) are inferior estimates. Hence, we also consider MPC2, where 
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Model-Predictive Control is used from the second stage of the conceptual algorithm. The 
first stage uses a predetermined sample size and number of solver iterations. We 
examine three choices for the first-stage policy resulting in the following three cases of 
MPC2; 

1. MPC2a. n, =450, = 100, and remaining stages use MPC. 

2. MPC2b. n, = 600, = 100, and remaining stages use MPC. 

3. MPC2c. = 900, = 100, and remaining stages use MPC. 

Choices of n, for these cases of MPC2 are determined by solving 6*"' = 0.01, 6*”' = 0.05, 

and 6*"' =0.1, where 6 is the conservative rate-of-convergence coefficient as discussed 
in Chapter 3, Section B, Subsection 2. 

As a basis for comparison, we consider the following heuristic policies: 

1. Fixed policy. Predetermine and n^. and keep fixed throughout each 
stage of CA. 

2. Additive policy. Predetermine and add a predetermined number to N^. 

at the beginning of each stage of CA (i.e., = 50 -i- A^.). 

3. Multiplicative policy. As in additive policy, predetermine and adjust 
sample size by a multiplicative factor at the beginning of each stage 
(i.e.,A,,,=1.2A,). 

We use an initial sample size Aj =10 for all heuristic policies, except for the 

fixed policy for which A^ =0.5'A* =5000 for all k. In order to use the PGM within 

CA, we must first find an initial feasible solution for both SCF and MCF to start the 
calculations. We do so by formulating and solving the following linear program in the 
case of SCF: 
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SCF-LP; min V C X 

i—4 ‘J V 

-D if / = 5 

s.t. Yj ^ij=' 0 iii^N\{s,t} 

KJaX^A ,:0J)^A ^ ^ 

0<Xy<M.j, \f i,j 

and the following linear program in the ease of MCF: 

MCF-LP; 

(iJ)^A p 

-D” if/ = 5 

s.t. H ^Ji~ H 0 if/e A^\{5,t}, Vp 

KJa)^a .(ij).A 

^x/’<M.., y i,j 

P 

0<x/’ , V i,j,p 

We implement our network-flow problems in Matlab Version 7.7.0 on a desktop 
computer running Windows XP with 3.73 GHz processor speed and 3.25 GB of RAM. 
SCF-LP and MCF-LP are solved to find an initial feasible solution for both SCF and 
MCF using the linear programming solver linprog in the optimization toolbox. 

For comparison studies, we record the computing time of CA using each of the 
MPC policies with the computing time for each of the other policies considered. We 
evaluate each of the heuristic policies with a different predetermined control on the 
number of solver iterations. Evaluations are run with for all k set at 5, 25, 50, 75 and 

100 iterations. In the additive approach, the sample size is increased by 100 at the 
beginning of every stage. For the multiplicative approach, two separate heuristics are 
considered. First, we evaluate the policy with an adjustment to sample size as 
A^^j=L5Aj. for all k and then increase the adjustment control on sample size to 

= 2A^ for all k. 
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The results summarized in Table 1 provide average eomputing times over 20 runs 
of the CA with standard deviations for the SCF problem. The first eolumn lists the 
individual polieies mentioned above for determining (A^, ). The seeond and third 

eolumns give the average and standard deviations, respeetively, of the total 
eomputational times to reaeh a near-optimal objeetive value in the SCF problem. 


Policy 

SCF Computational Times (sec.) 

avg st dev 

MFCl 

23.17 

3.60 

MFC2a,ni = 450, Ai = 100 

14.87 

2.67 

MFC2b,ni = 600, Ai = 100 

18.38 

2.48 

MFC2e,ni = 900, Ai = 100 

24.77 

3.15 

Fixed, n = 5 

443.30 

2.48 

Fixed, n =25 

638.12 

42.84 

Fixed, n = 50 

605.40 

50.24 

Fixed, n =15 

614.86 

38.27 

Fixed, n = 100 

631.84 

50.60 

Additive, n = 5 

398.07 

16.78 

Additive, n =25 

87.58 

7.78 

Additive, n = 50 

49.26 

6.76 

Additive, n =15 

41.54 

6.81 

Additive, n = 100 

40.81 

9.08 

Mult 1.5, n =5 

> 1100 


Mult 1.5, n = 25 

81.95 

13.54 

Mult 1.5, n = 50 

31.20 

6.31 

Mult 1.5, n = 75 

29.85 

7.85 

Mult 1.5, n = 100 

27.95 

5.93 

Mult 2.0, n = 5 

> 1100 


Mult 2.0, n =25 

> 1100 


Mult 2.0, n = 50 

77.06 

7.52 

Mult 2.0, n =15 

35.98 

7.34 

Mult 2.0, n = 100 

32.56 

8.29 


Table 1. Average and standard deviation of eomputing times of CA with Model-Predietive 
Control (MFC) polieies and heuristie polieies applied to SCF. 


In Table 1, the best heuristie poliey with respeet to eomputational time is the 


multiplieative poliey with a multiplieative faetor of 1.5 and n^=100 for all k. In 

eomparison, MFC I finds a near-optimal objeetive value nearly five seeonds faster and 
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does so with a 39.3% reduction in standard deviation of computational time over the 20 
independent runs. MPC2a improves further still, offering a reduction in computing time 
of nearly 47%. Additionally, the standard deviation between the independent runs drops 
55% as compared to the best heuristic policy. The computational times recorded do not 
reflect the time required to determine the Model-Predictive Control, i.e., to solve 
approximately S-SSCPj.(r^,r*,rg,r^). We elect to exclude this time because for large, 

real-world problems, computing times for the minimization calculations and checking 
stopping criterion are expected to be considerably larger than computing times for 
determining (A^,n^). 

Several policies considered for SCF return results that are costly regarding 
computing times. In those cases, we terminate the calculations after 1100 seconds and do 
not compute averages: see rows 12, 17, and 18 of Table 1. For the policies with times 
greater than 1100 seconds, the relatively small number of solver iterations is not 
sufficient to make substantial gains towards /*. In these cases, grows quite large 
and computing time suffers from the large sample size. For each of the problems, is 

limited to 400,000 to avoid exhausting computer memory, and in each of these cases, the 
sample size grows to this limit. 

The results summarized in Table 2 for MCF, provide average computing times 
over 20 runs of the CA, along with standard deviations. 
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Policy 

MCF Computational Times (sec.) 

avg std dev 

MPCl 

40.70 

7.21 

MPC2a, ni = 450, Ai = 100 

26.46 

4.22 

MPC2b,ni = 600, Ai = 100 

28.84 

4.35 

MPC2c,ni = 900,A^i = 100 

37.57 

2.95 

Fixed, n = 5 

> 1000 


Fixed, n =25 

891.65 

47.96 

Fixed, n = 50 

905.67 

61.53 

Fixed, n =15 

905.11 

58.92 

Fixed, n = 100 

891.89 

57.60 

Additive, n = 5 

612.40 

41.05 

Additive, n =25 

143.12 

14.81 

Additive, n = 50 

79.86 

9.91 

Additive, n =15 

66.74 

10.75 

Additive, n = 100 

54.86 

7.24 

Mult 1.5, n =5 

> 1000 


Mult 1.5, n = 25 

151.96 

33.66 

Mult 1.5, n = 50 

48.66 

8.86 

Mult 1.5, n = 75 

44.41 

7.51 

Mult 1.5, n = 100 

48.12 

12.71 

Mult 2.0, n = 5 

> 1000 


Mult 2.0, n =25 

> 1000 


Mult 2.0, n = 50 

142.55 

46.12 

Mult 2.0, n =15 

61.95 

16.52 

Mult 2.0, n = 100 

52.64 

11.67 


Table 2. Average and standard deviation of computing times of CA with Model-Predictive 
Control (MPC) policies and heuristic policies applied to MCF. 


As in SCF, a number of the policies make the sample size grow until it hits the 
limit of 400,000, thereby affecting overall computing times. In these cases we terminate 
the calculations after 1000 seconds and do not compute averages. 

The best heuristic policy for MCF is again a multiplicative policy. However in 
this larger problem, computational time is best when = 75 for all k: compare this to 

SCF, where computational time is best when =100 for all k. MPCl improves on this 
computational time by nearly four seconds and does so with essentially the same 

variability of computational time between runs as the best heuristic policy. 
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We see that modifying MPC in the first stage gives further eomputational savings. 
Poliey MPC2a shows an improvement of 40% in overall eomputing time, on average, and 
improves the standard deviation of eomputing time by almost 44% over the 20 
independent runs. These results indieate that while the MPC typieally provides a “good” 
poliey for seleeting {N^, «^), the estimates of parameters in S-SSCP^(r^,r*,r^,r^) for 
the first stage are rather poor and a heuristic policy may be better in that stage. 

To verily that the stopping test (3.5) does not cause premature termination of CA, 
we compute a lower bound on f* as described in Mak et. ah, (1999). Specifically, we 
run the PGM on the AP with N = 10000 until that algorithm stalls and record the last 
function value. This is an estimate of . We repeat this process 30 times. By the 
central limit theorem, the average of these function values is approximately normal and 
provides a lower bound /on /*. We find that in all 160 runs of the MPC policies, the 
probability that the last solution found is no worse than (1 + 0.01) / is essentially 1.0. 

Hence, the stopping test (3.5) is rather conservative, as zero unsatisfactory solutions is 
well within the 0.05 T60 = 8 expected when S = 0.95. 

MPCl solves each of the network-flow problems faster than any of the heuristic 
policies considered. Additional reductions in computing time are gained with MPC2 
where the first stage of the CA is primarily used to estimate parameters in 
S-SSCP^(r^,r*,r^,r^) fork>2. 
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V. CONCLUSIONS 


This thesis develops an efficient optimization algorithm for approximately solving 
stochastic nonlinear programming problems whose objective functions are sample- 
average approximations. We demonstrate improvement in computation times by 
approximately solving a discrete-time optimal-control problem to select a policy of well- 
balanced sample sizes and number of solver iterations for each stage of the algorithm. 
This policy, referred to as the Model-Predictive Control policy (MPC), is compared 
against alternative heuristic policies for selecting sample sizes and solver iterations. 
MPC approximately solves a single-commodity network-flow problem up to 17% faster, 
on average, than the best heuristic policy. Furthermore, the optimal-control problem 
provides a 40% reduction in standard deviation of computing times over a set of 
independent runs of the algorithm on identical problem instances. When we fix the 
number of solver iterations in the first stage and then proceed with MPC, we improve the 
computing time, on average, by nearly 47% and reduce the standard deviation between 
runs by more than one half. 

The application of the discrete-time optimal-control problem to a larger multi- 
commodity network-flow problem shows an 8.4% improvement, on average, in 
computational time over the best heuristic policy with essentially the same variation of 
overall computational time between the 20 independent runs. With the first-stage 
modification to Model-Predictive Control, we improve computing time by 40%, on 
average, compared to the same heuristic policy and reduce standard deviation between 
runs by 44%. 

The algorithm developed to solve nonlinear stochastic programs shows 
considerable promise and offers significant potential for further study. 
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